#!/usr/bin/python
"""
This program pick primers for a list of given regions.

Usage: dissertation_PrimerPicker.py genome.fasta regions.list
Note: Need to run in primer3 program folder:
/home/js/tools/primer3-2.3.2

Regions need to be in tab-delimited format like this:

region1	1000	2000
region2	6000	10000
region3	12000	14000

Currently, the values are hard-coded for long PCR primers to produce
product size of about 15kb.
Should check the results manually to pick best primers but eventually this script
should be able to pick the best one for given conditions (TODO).

"""

import sys
from Bio import SeqIO
from Bio import Motif
from Bio.Alphabet import IUPAC
from Bio.Seq import Seq

genome = SeqIO.read(sys.argv[1], "fasta")

primerfile = open(sys.argv[2], "rU")
primers = primerfile.readlines()

primerdict = {}
lefthit = {}
righthit = {}

for pair in primers:
    l = pair.split('\t')
    name = l[0]
    leftprimer = l[1]
    rightprimer = l[2].rstrip()
    left = Motif.Motif(alphabet=IUPAC.unambiguous_dna)
    left.add_instance(Seq(leftprimer, left.alphabet))
    right = Motif.Motif(alphabet=IUPAC.unambiguous_dna)
    tmpseq = Seq(rightprimer, right.alphabet).reverse_complement()
    right.add_instance(Seq(rightprimer, right.alphabet).reverse_complement())
    primerdict[leftprimer] = name
    primerdict[tmpseq] = name

    #searching the genome for primer-binding sites
    for lpos, lseq in left.search_instances(genome.seq):
        print name + "-F" + "\t" + str(lpos) + "->" + str(lpos+len(leftprimer)) + "\t" + lseq.tostring()
        lefthit[name] = ((lpos, lseq.tostring()))
    for rpos, rseq in right.search_instances(genome.seq):
        print name + "-R" + "\t" + str(rpos) + "->" + str(rpos+len(rightprimer)) + "\t" + rseq.tostring()
        righthit[name] = ((rpos, rseq.tostring()))

#getting left and right hit coordinates to calculate size of PCR product
for k, v in lefthit.iteritems():
    leftname = k
    leftpos = v[0]
    tmp = righthit[k]
    rightpos = tmp[0]
    productsize = int(rightpos) - int(leftpos) + 1
    print 'Region {0:4} starts at {1:7} and ends at {2:7} with product size of {3:5} bp'.format(leftname, leftpos, rightpos+len(tmp[1]), str(productsize))

primerfile.close()
